\documentclass[fleqn,english]{article}               % final
\usepackage[T1]{fontenc}
\usepackage{booktabs}			
\usepackage{multirow}
\usepackage{babel}
\usepackage{url}
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[margin=1.5in]{geometry}
\usepackage{multirow}
\usepackage[verbose]{placeins}

\title{Mathematical Programming \\ Programming exercise}

\author{
  Christian Hafner
  \thanks{christianhafner@gmail.com}
    \and
  Klemens Jahrmann
  \thanks{klemens.jahrmann@tuwien.ac.at}
  \and
  Kevin Streicher
  \thanks{kevin.streicher.at@gmail.com}
}

\begin{document}

\maketitle

\section{Notation}

Given is a weighted undirected graph $G'=(V,E',f')$ with $V=\{v_0,v_1,\ldots,v_{n-1}\}$ and $f':E' \rightarrow \mathbb{N}_0$. $v_0$ is a virtual node which is adjacent to all other nodes via edges of weight 0. The goal is to find a $k$-MST of $\{v_1,\ldots,v_{n-1}\}$, i.e.\ a subgraph of $G'$ which is a tree of $k$ nodes and has minimum weight.

We transform $G'$ into a directed graph $G = (V,E,f)$. Every edge $v_0w$ is replaced by a directed edge $(v_0,w)$. Every edge $vw$ not incident to $v_0$ is replaced by two directed edges $(v,w)$ and $(w,v)$. Every directed edge is assigned the weight of its corresponding undirected edge. The symbols $\Gamma^+$ and $\Gamma^-$ are used with respect to $G$.

The task of finding a $k$-MST in $G'$ is equivalent to finding an arborescence in $G$ by selecting a subset of $E$.

We based our solution for the CEC on the lecture slides and for the DCUT on the proposed formulation by Chimani, Kandyba and Ljubi\'c\footnote{\url{http://www.researchgate.net/publication/220982068_Obtaining_Optimal_k-Cardinality_Trees_Fast/file/79e4150eabae642f98.pdf}}. Also we used their \textbf{orientation} and \textbf{asymmetry} constraints, which improved the runtime down to our final results.

\section{Arborescence}
There are several rules which define our arborescence. These rules are used for both formulations.

\begin{enumerate}
\item We select exactly one of the outgoing edges of $v_0$, thereby choosing a unique root:
\begin{equation}
\sum_{v \in V \setminus \{v_0\} } e_{v_0 v} = 1.
\label{eq:c1}
\end{equation}

\item A node can have at most one incoming edge. This constraint forbids nodes with more than one parent. We denote $I_v$ as the sum of incoming edges for node $v$.
\begin{equation}
\forall v \in V: \sum_{u \in \Gamma^-(v)} e_{uv} \leq 1.
\label{eq:c2}
\end{equation}

\item A node can only have outgoing edges if it has one incoming edge. To reduce the Big-M as much as possible we use the maximum number of total edges as M, which is equal to $(k-1)$.
\begin{equation}
\forall v \in V: \sum_{u \in \Gamma^+(v)} e_{uv} \leq (k-1) \sum_{u \in \Gamma^-(v)} e_{uv}.
\label{eq:c3}
\end{equation}

\item The last constraint ensures that the tree has exactly $k$ nodes by using the property of trees $|V| = |E|+1$:
\begin{equation}
\sum_{(v,w) \in E} e_{vw} = k.
\label{eq:c4}
\end{equation}
This includes the $k-1$ edges of the $k$-MST and the virtual edge from $v_0$ to the root of the $k$-MST.

\item Each node with an incoming edge has to be selected:
\begin{equation}
\forall v \in V: z_v = \sum_{u \in \Gamma^-(v)} e_{uv}.
\label{eq:c5}
\end{equation}

\item The objective function is always the minimum sum of selected edge weights:
\begin{equation}
\min \sum_{\forall (u,v) \in E} w_{uv} \; e_{uv}.
\label{eq:c6}
\end{equation}

\item We added the strengthening constraints for orientability of arcs (\textbf{orientation-constraint}).
\[\forall i \in V,\forall\{i,j\} \in E: x_{ij}+x_{ji} \leq z_i\]

\item To remove symmetric solutions, we added the asymmetric-constraints. Lubi\'c et al. note, that these constraints might have a bad impact for bigger instances, for our largest provided instances it improved the runtime by an order of magnitude:
\[\forall i,j \in V, i<j : x_{0j} \leq 1 - z_{i}\]
\end{enumerate}

\section{Cycle Elimination Constraints}

In order to ensure acyclicity and connectedness, it suffices to forbid cycles. Connectedness follows from constraints~\eqref{eq:c3} and~\eqref{eq:c4}.

The exponential number of constraints that ensure acyclicity can be formulated as
\begin{equation}
	\forall C \subseteq E, |C| \geq 2, C \; \textrm{forms a cycle:} \sum_{(u,v) \in C} e_{uv} \leq |C|-1.
\end{equation}
These constraints are realized in a branch-and-cut manner, i.e. whenever an optimal solution satisfying the current constraints is found, we identify a violated cycle constraint and add it to the model, if such a violation exists.

This is implemented by considering the following transformation:
\begin{equation}
	\sum_{(u,v) \in C} e_{uv} \leq |C|-1 \iff \sum_{(u,v) \in C} (1-e_{uv}) \geq 1.
\end{equation}
Every arc $(u,v)$ with $e_{uv} = 1$, we look for the shortest path from $v$ to $u$. If this shortest path is 0, we can conclude that a cycle including the edge $(u,v)$ is part of our incumbent solution. This cycle is forbidden by adding the corresponding constraint. We add a multiple of $\varepsilon$ for numerical stability:
\begin{equation}
	\sum_{(u,v) \in C} e_{uv} \leq |C|(1+\varepsilon).
\end{equation}


\section{Directed cutset constraints}

In this formulation we ensure connectedness rather than acyclicity to arrive at a valid solution. This can be ensured by an exponential number of constraints that were formulated in  in the aforementioned work by Chimani et al.:
\begin{equation}
	\forall T \subseteq V \setminus \{v_0\}, \forall v \in T: \; \sum_{(s,t) \in \delta^-(T)} e_{st} \geq z_v,
\end{equation}
where $\delta^-(T)$ denote all arcs $(s,t)$ with $s \in V \setminus T$ and $t \in T$. These constraints ensure that there is a selected arc between any node set containing $v_0$ and its complementary node set.

At each separation step we set the capacities of the arcs according to their corresponding arc variables. Then we compute the maximum flow between $v_0$ and every other node $v$. When the maximum flow is lower than the value of $v$, we have detected a violation of connectedness and can derive the corresponding constraint from the cutset.

It was unexpected that this formulation yielded correct results for every instance without having to use $\varepsilon$ to ensure numerical stability.

\section{Results}
% Please remember to add \use{multirow} to your document preamble in order to suppor multirow cells

\begin{table}[h!]
\hspace{-55px}
    \begin{tabular}{lrrr|rrrr|rrrr}
    \toprule
          &       &       &       & \multicolumn{2}{c}{CEC} &       &       & \multicolumn{2}{c}{DCC} &       &  \\
    \midrule
    Instance & Nodes & k     & Optimal & CPU time[s] & Value & Nodes & Cuts  & CPU time[s] & Value & Nodes & Cuts \\
    g01   & 10    & 2     & 46    & 0,015 & 46    & 0     & 0     & 0,062 & 46    & 0     & 0 \\
          &       & 5     & 477   & 0,015 & 477   & 0     & 0     & 0,015 & 477   & 0     & 0 \\    \midrule
    g02   & 20    & 4     & 373   & 0,062 & 373   & 0     & 0     & 0,046 & 373   & 0     & 2 \\
          &       & 10    & 1390  & 0,062 & 1390  & 0     & 2     & 0,015 & 1390  & 0     & 7 \\    \midrule
    g03   & 50    & 10    & 725   & 0,015 & 725   & 0     & 2     & 0,062 & 725   & 0     & 19 \\
          &       & 25    & 3074  & 0,218 & 3074  & 188   & 26    & 0,171 & 3074  & 2     & 117 \\    \midrule
    g04   & 70    & 14    & 909   & 0,125 & 909   & 0     & 0     & 0,156 & 909   & 0     & 8 \\
          &       & 35    & 3292  & 0,078 & 3292  & 0     & 2     & 0,093 & 3292  & 0     & 5 \\  \midrule
    g05   & 100   & 20    & 1235  & 0,328 & 1235  & 5     & 7     & 0,171 & 1235  & 0     & 18 \\
          &       & 50    & 4898  & 0,203 & 4898  & 13    & 9     & 0,125 & 4898  & 0     & 84 \\\midrule
    g06   & 200   & 40    & 2068  & 7,515 & 2068  & 2044  & 103   & 2,156 & 2068  & 24    & 510 \\
          &       & 100   & 6705  & 3,343 & 6705  & 402   & 79    & 1,109 & 6705  & 0     & 454 \\\midrule
    g07   & 300   & 60    & 1335  & 3,578 & 1335  & 0     & 4     & 2,593 & 1335  & 0     & 13 \\
          &       & 150   & 4534  & 6,937 & 4534  & 13    & 8     & 2,671 & 4534  & 0     & 88 \\\midrule
    g08   & 400   & 80    & 1620  & 6,125 & 1620  & 0     & 2     & 5,687 & 1620  & 0     & 90 \\
          &       & 200   & 5787  & 8,687 & 5787  & 0     & 6     & 4,686 & 5878  & 0     & 97 \\
    \bottomrule
    \end{tabular}%
      \caption{Final results of our implemnetation for CEC and DCUT}
  \label{tab:addlabel}%
\end{table}%



\FloatBarrier
\section{Discussion}
For the small ($n<100$) instances both formulations are so fast, it is not really possible to make conclusions about them on the given instances, because a simple reordering of the constraints in cplex would have a bigger impact.

As two of us already implemented the MCF,SCF,MTZ formulations in Algorithmics, it is interesting to state that the CEC and DCC formulations are way faster. Exact comparisons cannot be made because the strengthening \textbf{orientation-} and \textbf{asymmetric-constraint} could have been implemented there as well and they already improved the CEC and DCC formulations a lot.

Considering the number of branch and bound nodes, both formulations need many fewer branches, which explains why they are so fast. The DCC formulation could solve 14 out of 16 instances without any branching at all and more violated constraints than for CEC were added. CEC could only solve 10 out of 16 instances in the root node and the highest number of branch and bound nodes of 2044 is approximately 85 times larger than the highest number of branch and bound nodes for DCC which was 24. This conclusion obviously can not be generalized as we would have needed more test instances like the instance libraries used in literature to make assumptions about this.

We did not implement neither nested nor backward cuts, although they have been reported to improve the running time a lot. We do not expect them to have a significant impact on those small instances as Chimani et al. also tried their formulation on instances ten times as large as ours.

\end{document}